Este documento contiene el código de R empleado para la elaboración del proyecto final de series de tiempo.

1. Importar y analizar datos

Los datos corresponden a la serie del índice accionario Standard and Poor’s 500, y se obtuvieron del símbolo GSPC de acuerdo a la información histórica disponible en https://finance.yahoo.com/quote/%5EGSPC. Se empleó el precio ajustado de cierre del índice del 3 de enero de 2020 al 20 de abril de 2021.

Cargar las librerías requeridas:

library("readxl")
library("xts")
library("tseries")
library("forecast")
library("rugarch")
library("vars")
library("PerformanceAnalytics")
library("DMwR")

Importar datosy revisar datos:

datos <- read_excel("gspc.xlsx")
summary(datos)
##       date                         close       
##  Min.   :2000-01-03 00:00:00   Min.   : 676.5  
##  1st Qu.:2005-05-03 06:00:00   1st Qu.:1165.9  
##  Median :2010-08-26 12:00:00   Median :1394.8  
##  Mean   :2010-08-27 12:04:34   Mean   :1684.8  
##  3rd Qu.:2015-12-21 18:00:00   3rd Qu.:2085.4  
##  Max.   :2021-04-20 00:00:00   Max.   :4185.5

Crear el objeto sp500 como una serie de tiempo (un objeto de la clase xts):

sp500 <- xts(x = datos$close, order.by = datos$date)

Verificar que la periodicidad de la serie sea correcta (se esperan datos diarios de 20 años):

periodicity(sp500)
## Daily periodicity from 2000-01-03 to 2021-04-20

Efectivamente, la serie de tiempo es de periodicidad diaria, con datos del 3 de enero de 2020 al 20 de abril de 2021.

Se grafica la serie de tiempo para una primera inspección visual:

plot(sp500)

Se prueba estacionariedad con la prueba aumentada de Dickey-Fuller:

adf.test(sp500, alternative = 's')
## 
##  Augmented Dickey-Fuller Test
## 
## data:  sp500
## Dickey-Fuller = -0.49676, Lag order = 17, p-value = 0.982
## alternative hypothesis: stationary

El valor \(p\) cercano a uno indica que no se puede rechazar la hipótesis nula de que la serie presenta raíz unitaria.

Se calcula una serie que contiene las primeras diferencias para verificar si el proceso es integrado de orden uno. Esta serie se denomina sp500_1:

sp500_1 <- log(sp500) - log(lag(sp500))
sp500_1 <- na.omit(sp500_1)

Verificar periodicidad de la serie de primeras diferencias:

periodicity(sp500_1)
## Daily periodicity from 2000-01-04 to 2021-04-20

Se tiene una observación menos debido al cálculo de primeras diferencias.

Se repite la prueba aumentada de Dickey-Fuller con la serie de primeras diferencias:

adf.test(sp500_1, alternative = 's')
## Warning in adf.test(sp500_1, alternative = "s"): p-value smaller than printed p-
## value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  sp500_1
## Dickey-Fuller = -17.641, Lag order = 17, p-value = 0.01
## alternative hypothesis: stationary

El valor \(p\) menor a 0.01 indica que se rechaza la hipótesis nula de que la serie sp500_1 tiene raíz unitaria. Como esta serie contiene las primeras diferencias de sp500, entonces esta última es una serie integrada de orden 1.

Gráfica de la serie de primeras diferencias:

plot(sp500_1)

Autocorrelación de la serie sp500:

sp500_acf <- acf(sp500)

Autocorrelación parcial de la serie sp500:

sp500_pacf <- pacf(sp500)

Autocorrelación de la serie sp500_1:

sp500_1_acf <- acf(sp500_1)

Autocorrelación parcial de la serie sp500_1:

sp500_1_pacf <- pacf(sp500_1)

Se realiza la prueba de Ljung-Box para deterninar si existe correlación serial en la serie estacionaria (sp500_1), para órdenes desde 1 hasta 30:

for (i in 1:30) {
  prueba = Box.test(sp500_1, lag = i, type = "Ljung")
  cat(i, "\t", round(prueba$statistic, 1), "\t", prueba$p.value, "\n")
}
## 1     68.8    1.110223e-16 
## 2     68.9    1.110223e-15 
## 3     70.1    3.996803e-15 
## 4     73.3    4.551914e-15 
## 5     74.3    1.276756e-14 
## 6     81.9    1.44329e-15 
## 7     87.8    3.330669e-16 
## 8     92.5    1.110223e-16 
## 9     100.3   0 
## 10    100.3   0 
## 11    101.5   1.110223e-16 
## 12    113.8   0 
## 13    115     0 
## 14    115.1   0 
## 15    135.7   0 
## 16    166.5   0 
## 17    166.7   0 
## 18    176.3   0 
## 19    176.3   0 
## 20    176.3   0 
## 21    178     0 
## 22    178     0 
## 23    178.2   0 
## 24    178.4   0 
## 25    178.4   0 
## 26    179.6   0 
## 27    184.9   0 
## 28    185.1   0 
## 29    185.6   0 
## 30    186.6   0

Los valores \(p\) iguales a cero indican que se rechaza la hipótesis nula de que las observaciones no tienen correlación serial para órdenes menores o iguales a 30.

2. Modelo ARMA

De acuerdo a la documentación del comando ARMA, la especificación del modelo ARMA(\(p\), \(q\)) en R es

\[ y_t = a_0 + a_1 y_{t-1} + \cdots + a_p y_{t-p} + b_1 e_{t-1} + \cdots + b_q e_{t-q} + e_t\]

Los comandos para los criterios de información en R sólo funcionan para estimaciones hechas mediante máxima verosimilitud (MLE), y la implementación del comando arma sólo es capaz de estimar usando mínimos cuadrados ordinarios (OLS). El comando que emplea MLE es Arima, por lo que es el que se utilizará para ajustar los modelos y obtener criterios de información. Dado que se empleará la serie sp500_1 (que tiene removida la raíz unitaria presente en la serie original), el orden de integración del modelo debe ser cero (ARIMA(\(p\), \(0\), \(q\))), lo cual equivale a un modelo ARMA(\(p\), \(q\)).

Como primer ejercicio se estima el modelo ARMA(1, 1):

arima_1_0_1 <- Arima(sp500_1, order = c(1, 0, 1))
print(arima_1_0_1)
## Series: sp500_1 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1   mean
##       0.0071  -0.1225  2e-04
## s.e.  0.1140   0.1132  2e-04
## 
## sigma^2 estimated as 0.0001546:  log likelihood=15902.42
## AIC=-31796.83   AICc=-31796.82   BIC=-31770.49

Para determinar el mejor modelo empleando el criterio de información de Akaike (AIC), se estman 100 modelos ARMA(\(p\), \(q\)) diferentes, donde los parámetros de orden \(p\) y \(q\) varían independientemente en el rango de 1 hasta 10. Para cada modelo se obtiene el criterio de Akaike, y éstos se almacenan en la matriz criterios_aic, donde el número de fila corresponde al orden \(p\) y el número de columna corresponde al orden \(q\):

criterios_aic <- matrix(nrow = 10, ncol = 10)

for (p in 1:10) {
  for (q in 1:10) {
    modelo <- Arima(sp500_1, order = c(p, 0, q))
    criterios_aic[p,q] <- modelo$aic
  }
}

La matriz criterios_aic resulta ser:

print(criterios_aic)
##            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
##  [1,] -31796.83 -31794.85 -31797.92 -31799.09 -31797.62 -31805.47 -31804.77
##  [2,] -31794.97 -31793.88 -31792.37 -31798.09 -31797.17 -31815.28 -31802.85
##  [3,] -31793.83 -31791.83 -31793.82 -31801.75 -31799.75 -31815.01 -31811.51
##  [4,] -31798.27 -31800.00 -31801.23 -31799.93 -31804.20 -31815.45 -31815.01
##  [5,] -31797.20 -31797.85 -31799.79 -31805.14 -31824.08 -31822.41 -31819.64
##  [6,] -31808.80 -31814.52 -31815.80 -31816.47 -31813.69 -31821.95 -31829.05
##  [7,] -31807.15 -31804.92 -31811.22 -31815.32 -31815.08 -31829.58 -31826.89
##  [8,] -31806.75 -31820.05 -31801.15 -31826.01 -31824.20 -31822.00 -31823.08
##  [9,] -31808.56 -31830.77 -31825.19 -31820.85 -31824.14 -31833.60 -31832.09
## [10,] -31806.61 -31804.60 -31830.82 -31828.03 -31817.85 -31828.08 -31829.97
##            [,8]      [,9]     [,10]
##  [1,] -31803.20 -31804.44 -31802.50
##  [2,] -31818.27 -31802.54 -31800.50
##  [3,] -31815.39 -31817.91 -31817.64
##  [4,] -31809.93 -31810.88 -31823.46
##  [5,] -31821.80 -31830.96 -31841.80
##  [6,] -31820.85 -31827.64 -31836.98
##  [7,] -31835.50 -31830.42 -31834.89
##  [8,] -31831.50 -31831.22 -31834.37
##  [9,] -31831.05 -31844.89 -31840.76
## [10,] -31836.79 -31843.16 -31841.25

El valor mínimo es:

print(min(criterios_aic))
## [1] -31844.89

Que corresponde a la siguiente entrada:

which(criterios_aic == min(criterios_aic), arr.ind = TRUE)
##      row col
## [1,]   9   9

Por lo tanto, el modelo óptimo es el ARMA(9, 9). Se ajusta dicho modelo a los datos:

arima_9_0_9 <- Arima(sp500_1, order = c(9, 0, 9))
print(arima_9_0_9)
## Series: sp500_1 
## ARIMA(9,0,9) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2     ar3     ar4      ar5      ar6      ar7      ar8
##       0.1031  -0.3081  0.5003  0.3950  -0.2157  -0.5844  -0.1778  -0.4054
## s.e.  0.0706   0.0387  0.0796  0.0872   0.0602   0.0674   0.0805   0.0611
##          ar9      ma1     ma2      ma3      ma4     ma5     ma6     ma7     ma8
##       0.3980  -0.2165  0.3167  -0.5186  -0.3737  0.2616  0.5165  0.1524  0.3641
## s.e.  0.0694   0.0690  0.0623   0.0937   0.1013  0.0786  0.0806  0.0754  0.0649
##           ma9   mean
##       -0.3697  2e-04
## s.e.   0.0747  1e-04
## 
## sigma^2 estimated as 0.0001528:  log likelihood=15942.44
## AIC=-31844.89   AICc=-31844.73   BIC=-31713.17
residuales <- checkresiduals(arima_9_0_9)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(9,0,9) with non-zero mean
## Q* = 21.2, df = 3, p-value = 9.567e-05
## 
## Model df: 19.   Total lags used: 22

Valor \(p\) de la prueba de Ljung-Box sobre los residuales del modelo ajustado:

print(residuales$p.value)
## [1] 9.567302e-05

Debido a las que presenta la serie con periodicidad diaria, se convierte a datos mensuales empleando el valor de cierre (del último día del mes) como valor mensual. Esta nueva serie se denomina sp500m:

sp500m <- to.monthly(sp500)
sp500m <- sp500m$sp500.Close

Gráfica de la serie mensual:

plot(sp500m)

Prueba de Dickey-Fuller aumentada:

adf.test(sp500m, alternative = 's')
## Warning in adf.test(sp500m, alternative = "s"): p-value greater than printed p-
## value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  sp500m
## Dickey-Fuller = 0.24955, Lag order = 6, p-value = 0.99
## alternative hypothesis: stationary

No se puede rechazar la hipótesis nula de que la serie tiene raíz unitaria. Se obtiene la serie sp500m_1 que contiene las primeras diferencias:

sp500m_1 <- log(sp500m) - log(lag(sp500m))
sp500m_1 <- na.omit(sp500m_1)

Gráfica de sp500m_1:

plot(sp500m_1)

Prueba de DIckey-Fuller aumentada para sp500m_1:

adf.test(sp500m_1)
## Warning in adf.test(sp500m_1): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  sp500m_1
## Dickey-Fuller = -5.6307, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

Se rechaza la hipótesis nula de que la serie sp500m_1 tiene raíz unitaria.

Función de autocorrelación de sp500m_1:

sp500m_1_acf <- acf(sp500m_1)

Función de autocorrelación parcial de sp500m_1:

sp500m_1_pacf <- pacf(sp500m_1)

Prueba de Ljung-Box para determinar si existe autocorrelación serial:

for (i in 1:30) {
  prueba = Box.test(sp500m_1, lag = i, type = "Ljung")
  cat(i, "\t", round(prueba$statistic, 1), "\t", prueba$p.value, "\n")
}
## 1     1.5     0.2157794 
## 2     2.8     0.2453356 
## 3     4.4     0.2173632 
## 4     5.7     0.223933 
## 5     5.9     0.3124418 
## 6     7.5     0.2761677 
## 7     9.1     0.2486254 
## 8     9.6     0.2921704 
## 9     11.5    0.2440293 
## 10    11.6    0.3143923 
## 11    11.6    0.3937488 
## 12    11.7    0.4724017 
## 13    12.5    0.4898004 
## 14    13.4    0.4972188 
## 15    15.8    0.3935624 
## 16    16.9    0.3952379 
## 17    17      0.4555091 
## 18    17.6    0.484824 
## 19    21.3    0.319761 
## 20    23      0.287934 
## 21    23.2    0.3346201 
## 22    23.7    0.3619937 
## 23    23.7    0.4187952 
## 24    24.1    0.4574857 
## 25    26.3    0.3933144 
## 26    28.8    0.3217141 
## 27    29.2    0.3516093 
## 28    29.8    0.3710967 
## 29    30.4    0.392367 
## 30    30.4    0.4432944

Se repite la búsqueda del modelo ARMA óptimo para la nueva serie sp500m_1, esta vez variando los parámetros \(p\) y \(q\) entre 1 y 10:

criterios_aic_m <- matrix(nrow = 10, ncol = 10)

for (p in 1:10) {
  for (q in 1:10) {
    modelo <- Arima(sp500m_1, order = c(p, 0, q), method = 'ML')
    criterios_aic_m[p,q] <- modelo$aic
  }
}

La matriz criterios_aic_m resulta ser:

print(criterios_aic_m)
##            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
##  [1,] -867.2541 -865.3633 -864.7100 -862.7104 -861.8807 -860.7557 -859.5645
##  [2,] -865.3693 -863.3677 -862.7102 -864.6792 -866.0731 -864.1559 -861.7982
##  [3,] -864.6408 -866.7266 -860.7138 -862.9615 -864.1541 -862.1545 -860.1897
##  [4,] -862.7110 -864.7479 -863.3447 -866.0930 -864.2054 -860.1826 -863.0429
##  [5,] -862.2362 -865.7666 -865.8954 -864.2013 -862.1549 -858.2120 -857.1438
##  [6,] -861.5606 -863.7697 -861.8072 -859.7982 -858.1328 -856.1995 -860.2383
##  [7,] -860.7160 -861.5805 -859.8537 -860.7448 -860.7414 -859.1952 -862.6002
##  [8,] -859.7129 -860.6694 -858.0276 -858.8732 -858.1847 -862.6773 -855.7247
##  [9,] -858.6157 -856.6158 -854.6661 -859.8558 -857.7075 -860.8289 -858.8275
## [10,] -856.6165 -854.6157 -856.7950 -855.3184 -854.8144 -858.8293 -856.2785
##            [,8]      [,9]     [,10]
##  [1,] -859.6983 -859.3776 -857.9054
##  [2,] -860.6344 -858.6660 -864.0275
##  [3,] -858.2261 -856.6941 -861.3563
##  [4,] -860.0398 -854.6932 -854.4539
##  [5,] -858.3866 -857.8570 -855.7522
##  [6,] -860.4172 -860.8915 -858.9088
##  [7,] -857.2897 -858.8627 -856.9329
##  [8,] -857.2197 -857.0902 -855.2309
##  [9,] -857.1309 -857.0570 -856.0372
## [10,] -854.7100 -855.8460 -859.5798

El valor mínimo es:

print(min(criterios_aic_m))
## [1] -867.2541

Que corresponde a la siguiente entrada:

which(criterios_aic_m == min(criterios_aic_m), arr.ind = TRUE)
##      row col
## [1,]   1   1

Por lo tanto, el modelo óptimo es el ARMA(1, 1). Se ajusta dicho modelo a los datos:

arima_m_1_0_1 <- Arima(sp500m_1, order = c(1, 0, 1))
print(arima_m_1_0_1)
## Series: sp500m_1 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##           ar1     ma1    mean
##       -0.6531  0.7552  0.0043
## s.e.   0.2674  0.2323  0.0029
## 
## sigma^2 estimated as 0.001914:  log likelihood=437.63
## AIC=-867.25   AICc=-867.09   BIC=-853.09

Residuales:

residuales_m <- checkresiduals(arima_m_1_0_1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 22.279, df = 21, p-value = 0.3836
## 
## Model df: 3.   Total lags used: 24

Valor \(p\) de la prueba de Ljung-Box sobre los residuales del modelo ajustado:

print(residuales_m$p.value)
## [1] 0.3835919

Por tanto los residuales parecen ser ruido blanco.

Pronóstico dentro de la muestra:

training_set <- sp500m_1['/2018-12-31']
evaluation_set <- sp500m_1['2019-01-01/']

training_model <- Arima(training_set, order = c(1, 0, 1))

pronostico <- xts(predict(training_model, length(evaluation_set))$pred, order.by = index(evaluation_set))
pronostico_in <- merge(pronostico, evaluation_set, join='inner')

plot.xts(pronostico_in, screens = factor(1, 1),  legend.loc = 'bottomright')

Medidas de error del pronóstico dentro de la muestra:

medidas_error = regr.eval(coredata(evaluation_set), coredata(predict(training_model, length(evaluation_set))$pred))
print(round(medidas_error, 4))
##    mae    mse   rmse   mape 
## 0.0449 0.0031 0.0561 1.0233

Pronóstico fuera de la muestra:

pronostico <- xts(predict(training_model, 12)$pred, order.by = seq(as.Date("2021-06-01"),length=12,by="months"))
pronostico_out <- merge(pronostico, sp500_1, join = 'outer', fill = 0)['2021-01-01/']

plot.xts(pronostico_out, screens = factor(1, 1),  legend.loc = 'bottomright')

Modelo de volatilidad GARCH

Se estiman 25 modelos GARCH diferentes variando los parámetros \(p\) y \(q\) entre 1 y 5. Para cada modelo se calcula el criterio de información bayesiano (BIC) y éstos se almacenan en la matriz criterios_bic, donde el número de fila corresponde al orden de \(p\) y el número de columna corresponde al orden de \(q\):

criterios_bic = matrix(nrow = 5, ncol= 5)

for (p in 1:5) {
  for (q in 1:5) {
    # Especificación del modelo:
    garch_spec <- ugarchspec(
      variance.model = list(
        model = "sGARCH",
        garchOrder = c(p, q)),
      mean.model = list(
        armaOrder = c(p, q)),
      distribution.model = "norm")
    # Ajuste del modelo:
    garch_model <- ugarchfit(spec = garch_spec, data = sp500_1)
    # Cálculo del BIC:
    criterios_bic[p, q] = infocriteria(garch_model)[2]
  }
}

Matriz de criterios BIC:

criterios_bic
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] -6.449218 -6.446946 -6.443827 -6.440720 -6.437653
## [2,] -6.449317 -6.446842 -6.443833 -6.440217 -6.437456
## [3,] -6.445972 -6.443891 -6.440953 -6.439616 -6.434335
## [4,] -6.442883 -6.440892 -6.437719 -6.435750 -6.432870
## [5,] -6.439802 -6.437584 -6.438893 -6.432912 -6.430124

El valor mínimo es:

print(min(criterios_bic))
## [1] -6.449317

Que corresponde a la siguiente entrada:

which(criterios_bic == min(criterios_bic), arr.ind = TRUE)
##      row col
## [1,]   2   1

Por lo que el modelo óptimo es el GARCH(2, 1). Este modelo se ajusta en el objeto garch_2_1:

garch_spec <- ugarchspec(
  variance.model = list(
    model = "sGARCH",
    garchOrder = c(2, 1)),
  mean.model = list(
    armaOrder = c(2, 1)),
  distribution.model = "norm")

garch_2_1 <- ugarchfit(spec = garch_spec, data = sp500_1)

Resultado del ajuste:

show(garch_2_1)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(2,1)
## Mean Model   : ARFIMA(2,0,1)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.000621    0.000075    8.3143 0.000000
## ar1     0.891415    0.014314   62.2750 0.000000
## ar2     0.043317    0.014708    2.9452 0.003228
## ma1    -0.956180    0.001249 -765.3050 0.000000
## omega   0.000003    0.000001    3.3912 0.000696
## alpha1  0.081623    0.013697    5.9591 0.000000
## alpha2  0.064236    0.016971    3.7849 0.000154
## beta1   0.834458    0.012874   64.8185 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.000621    0.000076     8.1919 0.000000
## ar1     0.891415    0.013650    65.3032 0.000000
## ar2     0.043317    0.014219     3.0464 0.002316
## ma1    -0.956180    0.000397 -2408.4911 0.000000
## omega   0.000003    0.000004     0.6531 0.513690
## alpha1  0.081623    0.024367     3.3498 0.000809
## alpha2  0.064236    0.035469     1.8111 0.070133
## beta1   0.834458    0.044652    18.6880 0.000000
## 
## LogLikelihood : 17308.84 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.4592
## Bayes        -6.4493
## Shibata      -6.4592
## Hannan-Quinn -6.4557
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       1.072  0.3005
## Lag[2*(p+q)+(p+q)-1][8]      4.051  0.7620
## Lag[4*(p+q)+(p+q)-1][14]     7.657  0.4336
## d.o.f=3
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       1.211  0.2712
## Lag[2*(p+q)+(p+q)-1][8]      3.595  0.5770
## Lag[4*(p+q)+(p+q)-1][14]     7.396  0.4557
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]    0.5986 0.500 2.000  0.4391
## ARCH Lag[6]    2.0371 1.461 1.711  0.4826
## ARCH Lag[8]    2.8568 2.368 1.583  0.5706
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  12.6708
## Individual Statistics:              
## mu     0.70986
## ar1    0.12406
## ar2    0.09757
## ma1    0.18289
## omega  0.71000
## alpha1 0.19492
## alpha2 0.56190
## beta1  0.78566
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.89 2.11 2.59
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value      prob sig
## Sign Bias           3.2301 1.245e-03 ***
## Negative Sign Bias  0.2194 8.264e-01    
## Positive Sign Bias  2.3019 2.138e-02  **
## Joint Effect       39.5560 1.323e-08 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     196.6    1.638e-31
## 2    30     226.0    2.143e-32
## 3    40     259.8    2.035e-34
## 4    50     298.9    1.452e-37
## 
## 
## Elapsed time : 0.95836

Se observa que todos los coeficientes son significativos.

Se analiza qué sucede al aumentar el número de parámetros a \(p=3\) y \(q=2\):

garch_spec <- ugarchspec(
  variance.model = list(
    model = "sGARCH",
    garchOrder = c(3, 2)),
  mean.model = list(
    armaOrder = c(3, 2)),
  distribution.model = "norm")

garch_3_2 <- ugarchfit(spec = garch_spec, data = sp500_1)

show(garch_3_2)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(3,2)
## Mean Model   : ARFIMA(3,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.000621    0.000072    8.56391 0.000000
## ar1     0.594759    0.014221   41.82339 0.000000
## ar2     0.304970    0.016806   18.14609 0.000000
## ar3     0.016698    0.014684    1.13714 0.255482
## ma1    -0.659325    0.001306 -504.96046 0.000000
## ma2    -0.285140    0.003228  -88.34446 0.000000
## omega   0.000005    0.000001    3.99799 0.000064
## alpha1  0.079003    0.013758    5.74248 0.000000
## alpha2  0.143467    0.013550   10.58773 0.000000
## alpha3  0.025565    0.017609    1.45182 0.146552
## beta1   0.069598    0.092797    0.75001 0.453251
## beta2   0.649545    0.087091    7.45821 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.000621    0.000075    8.29435 0.000000
## ar1     0.594759    0.013379   44.45413 0.000000
## ar2     0.304970    0.017513   17.41380 0.000000
## ar3     0.016698    0.014742    1.13269 0.257345
## ma1    -0.659325    0.001295 -509.27693 0.000000
## ma2    -0.285140    0.000836 -341.14302 0.000000
## omega   0.000005    0.000005    1.00796 0.313473
## alpha1  0.079003    0.025352    3.11620 0.001832
## alpha2  0.143467    0.018735    7.65778 0.000000
## alpha3  0.025565    0.020267    1.26141 0.207160
## beta1   0.069598    0.138611    0.50211 0.615590
## beta2   0.649545    0.164953    3.93775 0.000082
## 
## LogLikelihood : 17311.48 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.4586
## Bayes        -6.4439
## Shibata      -6.4587
## Hannan-Quinn -6.4535
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       1.084  0.2977
## Lag[2*(p+q)+(p+q)-1][14]     7.572  0.4419
## Lag[4*(p+q)+(p+q)-1][24]    15.389  0.1388
## d.o.f=5
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       1.638  0.2006
## Lag[2*(p+q)+(p+q)-1][14]     5.989  0.6404
## Lag[4*(p+q)+(p+q)-1][24]     9.424  0.7595
## d.o.f=5
## 
## Weighted ARCH LM Tests
## ------------------------------------
##              Statistic Shape Scale P-Value
## ARCH Lag[6]      1.356 0.500 2.000  0.2443
## ARCH Lag[8]      1.851 1.480 1.774  0.5471
## ARCH Lag[10]     5.696 2.424 1.650  0.2136
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  6.2489
## Individual Statistics:              
## mu     0.66619
## ar1    0.13507
## ar2    0.06704
## ar3    0.21539
## ma1    0.19551
## ma2    0.13707
## omega  0.17338
## alpha1 0.14043
## alpha2 0.42789
## alpha3 0.74637
## beta1  0.77642
## beta2  0.82597
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.69 2.96 3.51
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value      prob sig
## Sign Bias           3.2132 1.320e-03 ***
## Negative Sign Bias  0.3199 7.491e-01    
## Positive Sign Bias  2.2909 2.201e-02  **
## Joint Effect       39.9379 1.098e-08 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     197.3    1.184e-31
## 2    30     230.7    2.673e-33
## 3    40     258.5    3.601e-34
## 4    50     286.6    2.598e-35
## 
## 
## Elapsed time : 1.314727

Normalidad de los residuales:

plot(garch_2_1, which = 9)

Los residuales no parecen seguir una distribución normal.

Autocorrelación de los residuales:

plot(garch_2_1, which = 10)

Gráfico de la volatilidad con 2 desviaciones estándar:

plot(garch_2_1, which = 1)

El valor en riesgo de la serie es (para \(p\) igual a 0.90, 0.95 y 0.99):

for (prob in c(0.90, 0.95, 0.99)) {
  print(c(prob, VaR(sp500_1, p = prob, method = 'historical')))
}
## [1]  0.90000000 -0.01301216
## [1]  0.95000000 -0.01913197
## [1]  0.99000000 -0.03520641

Modelo de vectores autorregresivos

Serie del retorno de largo plazo del Tesoro de EUA (compuesto a 10 años): https://www.treasury.gov/resource-center/data-chart-center/interest-rates/pages/TextView.aspx?data=longtermrateAll

Importar datos y revisar datos:

datos_ltc <- read_excel("lt_composite.xlsx")
summary(datos_ltc)
##       date                      lt_composite  
##  Min.   :2000-01-03 00:00:00   Min.   :0.970  
##  1st Qu.:2005-04-28 18:00:00   1st Qu.:2.740  
##  Median :2010-08-16 12:00:00   Median :3.810  
##  Mean   :2010-08-17 03:47:58   Mean   :3.798  
##  3rd Qu.:2015-12-07 06:00:00   3rd Qu.:4.820  
##  Max.   :2021-03-31 00:00:00   Max.   :6.910

Crear el objeto ltc como una serie de tiempo (un objeto de la clase xts):

ltc <- xts(x = datos_ltc$lt_composite, order.by = datos_ltc$date)

Verificar que la periodicidad de la serie sea correcta (se esperan datos diarios de 20 años):

periodicity(ltc)
## Daily periodicity from 2000-01-03 to 2021-03-31

Gráfica de la serie de tiempo:

plot(ltc)

Se prueba estacionariedad con la prueba ADF:

adf.test(ltc, alternative = 's')
## Warning in adf.test(ltc, alternative = "s"): p-value smaller than printed p-
## value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ltc
## Dickey-Fuller = -4.3485, Lag order = 17, p-value = 0.01
## alternative hypothesis: stationary

El valor \(p\) de 0.01 indica que se rechaza la hipótesis nula de que la serie presenta raíz unitaria. Aun así, para que esta serie sea consistente con la sp500_1, se laculan las primeras diferencias en una nueva serie ltc_1:

ltc_1 <- ltc - lag(ltc)
ltc_1 <- na.omit(ltc_1)

Verificar periodicidad de la serie de primeras diferencias:

periodicity(ltc_1)
## Daily periodicity from 2000-01-04 to 2021-03-31

Se tiene una observación menos debido al cálculo de primeras diferencias.

Se repite la prueba aumentada de Dickey-Fuller con la serie de primeras diferencias:

adf.test(ltc_1)
## Warning in adf.test(ltc_1): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ltc_1
## Dickey-Fuller = -15.967, Lag order = 17, p-value = 0.01
## alternative hypothesis: stationary

La serie conserva su estacionariedad original.

Gráfico de la serie de primeras diferencias:

plot(ltc_1)

Se encuentra el inicio y el final comunes de ambas series:

inicio <- max(min(index(sp500_1)), min(index(ltc_1)))
final <- min(max(index(sp500_1)), max(index(ltc_1)))

Estas fechas son, respectivamente:

print(c(inicio, final))
## [1] "2000-01-03 18:00:00 CST" "2021-03-30 18:00:00 CST"
res <- merge(sp500, ltc, join='inner')
VARselect(res)
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##     10     10     10     10 
## 
## $criteria
##                1         2         3         4         5         6         7
## AIC(n) 0.1771518 0.1528343 0.1477548 0.1478244 0.1460910 0.1473937 0.1343622
## HQ(n)  0.1797549 0.1571728 0.1538286 0.1556336 0.1556355 0.1586736 0.1473775
## SC(n)  0.1846005 0.1652488 0.1651351 0.1701705 0.1734028 0.1796713 0.1716056
## FPE(n) 1.1938123 1.1651320 1.1592287 1.1593094 1.1573015 1.1588101 1.1438071
##                8         9        10
## AIC(n) 0.1229181 0.1187283 0.1129452
## HQ(n)  0.1376687 0.1352143 0.1311666
## SC(n)  0.1651273 0.1659033 0.1650860
## FPE(n) 1.1307918 1.1260640 1.1195707
var_10 <- VAR(res, p = 10, type = 'const')
summary(var_10)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: sp500, ltc 
## Deterministic variables: const 
## Sample size: 5296 
## Log Likelihood: -15286.476 
## Roots of the characteristic polynomial:
##     1 0.9975 0.8676 0.8676 0.7582 0.7582 0.755 0.755 0.7129 0.6897 0.6897 0.6869 0.6869 0.6728 0.6728 0.6495 0.6495 0.6369 0.6369 0.01075
## Call:
## VAR(y = res, p = 10, type = "const")
## 
## 
## Estimation results for equation sp500: 
## ====================================== 
## sp500 = sp500.l1 + ltc.l1 + sp500.l2 + ltc.l2 + sp500.l3 + ltc.l3 + sp500.l4 + ltc.l4 + sp500.l5 + ltc.l5 + sp500.l6 + ltc.l6 + sp500.l7 + ltc.l7 + sp500.l8 + ltc.l8 + sp500.l9 + ltc.l9 + sp500.l10 + ltc.l10 + const 
## 
##            Estimate Std. Error t value Pr(>|t|)    
## sp500.l1    0.87675    0.01448  60.543  < 2e-16 ***
## ltc.l1     -1.87876    5.44382  -0.345 0.730020    
## sp500.l2    0.17905    0.01933   9.261  < 2e-16 ***
## ltc.l2     -9.25540    7.72063  -1.199 0.230664    
## sp500.l3   -0.02093    0.01935  -1.081 0.279546    
## ltc.l3      7.64278    7.71843   0.990 0.322122    
## sp500.l4   -0.08706    0.01915  -4.545 5.61e-06 ***
## ltc.l4      3.69813    7.71172   0.480 0.631569    
## sp500.l5    0.04667    0.01914   2.439 0.014778 *  
## ltc.l5     -4.34354    7.70504  -0.564 0.572964    
## sp500.l6   -0.09209    0.01915  -4.810 1.55e-06 ***
## ltc.l6     16.75382    7.70438   2.175 0.029706 *  
## sp500.l7    0.19857    0.01916  10.366  < 2e-16 ***
## ltc.l7    -28.94440    7.70762  -3.755 0.000175 ***
## sp500.l8   -0.15877    0.01939  -8.190 3.24e-16 ***
## ltc.l8     18.66746    7.71211   2.421 0.015531 *  
## sp500.l9    0.14443    0.01936   7.460 1.01e-13 ***
## ltc.l9     -6.24425    7.71494  -0.809 0.418338    
## sp500.l10  -0.08713    0.01450  -6.011 1.97e-09 ***
## ltc.l10     2.99994    5.43673   0.552 0.581115    
## const       4.75519    1.98274   2.398 0.016506 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 20.41 on 5275 degrees of freedom
## Multiple R-Squared: 0.9992,  Adjusted R-squared: 0.9992 
## F-statistic: 3.208e+05 on 20 and 5275 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation ltc: 
## ==================================== 
## ltc = sp500.l1 + ltc.l1 + sp500.l2 + ltc.l2 + sp500.l3 + ltc.l3 + sp500.l4 + ltc.l4 + sp500.l5 + ltc.l5 + sp500.l6 + ltc.l6 + sp500.l7 + ltc.l7 + sp500.l8 + ltc.l8 + sp500.l9 + ltc.l9 + sp500.l10 + ltc.l10 + const 
## 
##             Estimate Std. Error t value Pr(>|t|)    
## sp500.l1  -1.181e-04  3.864e-05  -3.057  0.00225 ** 
## ltc.l1     1.005e+00  1.453e-02  69.183  < 2e-16 ***
## sp500.l2   1.292e-04  5.159e-05   2.504  0.01230 *  
## ltc.l2    -5.484e-02  2.060e-02  -2.662  0.00779 ** 
## sp500.l3  -4.919e-05  5.164e-05  -0.952  0.34089    
## ltc.l3     4.448e-02  2.060e-02   2.160  0.03085 *  
## sp500.l4  -4.633e-05  5.111e-05  -0.906  0.36471    
## ltc.l4     7.148e-04  2.058e-02   0.035  0.97229    
## sp500.l5   8.841e-05  5.107e-05   1.731  0.08347 .  
## ltc.l5    -1.864e-03  2.056e-02  -0.091  0.92778    
## sp500.l6   3.611e-05  5.109e-05   0.707  0.47967    
## ltc.l6     5.426e-03  2.056e-02   0.264  0.79186    
## sp500.l7  -2.949e-05  5.112e-05  -0.577  0.56408    
## ltc.l7     1.737e-03  2.057e-02   0.084  0.93271    
## sp500.l8  -1.946e-05  5.173e-05  -0.376  0.70680    
## ltc.l8    -3.623e-02  2.058e-02  -1.761  0.07835 .  
## sp500.l9   1.025e-04  5.166e-05   1.984  0.04730 *  
## ltc.l9     3.062e-02  2.059e-02   1.487  0.13703    
## sp500.l10 -9.559e-05  3.868e-05  -2.471  0.01349 *  
## ltc.l10    2.915e-03  1.451e-02   0.201  0.84076    
## const      1.021e-02  5.291e-03   1.931  0.05359 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.05447 on 5275 degrees of freedom
## Multiple R-Squared: 0.9982,  Adjusted R-squared: 0.9982 
## F-statistic: 1.462e+05 on 20 and 5275 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##          sp500      ltc
## sp500 416.6996 0.354484
## ltc     0.3545 0.002967
## 
## Correlation matrix of residuals:
##        sp500    ltc
## sp500 1.0000 0.3188
## ltc   0.3188 1.0000
plot(var_10)

arch(var_10)
## Warning: Function 'arch' is deprecated; use 'arch.test' instead.
## See help("vars-deprecated") and help("arch-deprecated") for more information.
## 
##  ARCH (multivariate)
## 
## data:  Residuals of VAR object x
## Chi-squared = 3086.2, df = 45, p-value < 2.2e-16
ir <- irf(var_10, seed = 1)
plot(ir)

ir$runs
## [1] 100
ir$irf
## $sp500
##          sp500        ltc
##  [1,] 20.41322 0.01736539
##  [2,] 17.86466 0.01504078
##  [3,] 19.12880 0.01469069
##  [4,] 19.50848 0.01375604
##  [5,] 18.39539 0.01204305
##  [6,] 18.56075 0.01236441
##  [7,] 16.69747 0.01310932
##  [8,] 18.86988 0.01347808
##  [9,] 17.06684 0.01249074
## [10,] 19.27205 0.01472273
## [11,] 18.90516 0.01422970
## 
## $ltc
##             sp500        ltc
##  [1,]  0.00000000 0.05162921
##  [2,] -0.09699875 0.05188628
##  [3,] -0.66037459 0.04932462
##  [4,] -0.77465724 0.04908652
##  [5,] -0.75664287 0.04898153
##  [6,] -0.98155554 0.04869442
##  [7,] -0.27452383 0.04873090
##  [8,] -1.18072949 0.04871283
##  [9,] -0.90665051 0.04697710
## [10,] -1.15325776 0.04664313
## [11,] -1.19386703 0.04671415